Possible explanation for star-crushing effect in binary neutron star simulations 
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A possible explanation is suggested for the controversial star-crushing effect seen in numerical 
simulations of inspiraling neutron star binaries by Wilson, Mathews and Marronetti (WMM). An 
apparently incorrect definition of momentum density in the momentum constraint equation used 
by WMM gives rise to a post-l-Newtonian error in the approximation scheme. We show by means 
of an analytic, post-l-Newtonian calculation that this error causes an increase of the stars' central 
densities which is of the order of several percent when the stars are separated by a few stellar radii, 
in agreement with what is seen in the simulations. 
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A controversial issue in the astrophysics community re- 
cently has been the claim by Wilson, Mathews and Mar- 
ronetti (WMM), based on numerical simulations, that 
inspiraling binary neutron stars are subject to a general- 
relativistic crushing force that cause them to individually 
collapse to black holes before they merge Such a 

crushing force, if it existed, would have profound impli- 
cations for current efforts to detect gravitational waves 
from such systems with LIGO, VIRGO and other ground 
based detectors. The WMM claim has been disputed by 
several researchers utilizing a variety of approximate an- 
alytical and numerical techniques , and recent inde- 
pendent numerical simulations using the same approxi- 
mation scheme as WMM show no crushing effect [§). In 
this paper we suggest an explanation for the star-crushing 
effect perceived by WMM. 

We start with the standard ADM equations. The met- 
ric is 



ds 2 



-(a 2 - fr(3 l )dt 2 + 2P,dx l dt + -,,f/.r', /.,-'. (1) 



so that the lapse function is a and the shift vector is (3 l 
The extrinsic curvature Kij is given by 



7y = -2aK i:j + A/3j + Djfc, 



(2) 



where A is the derivative operator associated with 7^ 
and dots denote derivatives with respect to t. The Hamil- 
tonian constraint is 



(3) 



/,' /\,,/\'' • K 2 = 16irp H , 



(3) 



where ^R is the Ricci scalar of 7^-, K = ^ K^, n a is the 
normal to the t = const surface given by np = ~a{dt)p 1 
and ph = T a pn a rfi . The momentum constraint is 



Di(K ij -7^'iT) = 8tt5 3 , 



(4) 



where S a = -h^T^n 1 and h afi = g aSi + n a n is the 
projection tensor. [Here Greek indices run over (0, 1, 2, 3) 
and Roman indices over (1, 2, 3).] Finally the trace of the 
space-space part of Einstein's equation is 



47raS = K + AAa - l D t K 



^R + ZK^K^+K 2 



(5) 



where S = h af3 T a P. 

The main elements of the WMM approximation 
scheme are as follows fl|-|| : (i) They use the standard 
perfect fluid equations to evolve the fluid in the back- 
ground metric (uh. The stress-energy tensor is 



T a p = (p + p)u a Uf} + pg a p, 



(6) 



where u a is the 4- velocity, p is the pressure and p is the 
energy density. The equations of motion are V a T a ^ = 
and V Q (nii Q ) = 0, where n is the baryon number density, 
(ii) They work in a co-rotating coordinate system of the 
form (|lj), so that the large r boundary condition on the 
shift vector is (3i(x^) — eijk^x k + 0(r°), where fP is 
the orbital angular velocity, (iii) They impose that the 
spatial metric 7^- be conformally flat, 7y = p 4 jij, where 
7y is flat and time independent. By decomposing the 
extrinsic curvature as Kij = Aij + Kjij/3 where Ay- 
is traceless, and combining with @ and the conformal 
flatness condition one gets 



K = =«£ + -A/3* 

(pa a 



and 



A; 



2a 



D t p j +D J (3 l --{D k f3 k ) llJ 



(7) 



(8) 



Using the relations (0) and (||), the Hamiltonian con- 
straint (Q), the momentum constraint (Q) and the dy- 
namical equation (||) can be written schematically as 



<p = F 1 [ip,0 i ,a], 
= F^,a }i p,cp], 
<p = F 3 [(p,(p,a,(3 k ], 



(9) 
(10) 
(11) 



for some functionals Fi, F^ and A- (iv) They use 
a quasi-equilibrium approximation scheme which means 
that they substitute (p = (p = into Eqs. (p[)-(pd|). (v) 
They substitute the maximal slicing condition K — into 
the resulting equations. This yields a system of equations 
in which one can solve for a, ip, and (3 l at each instant 
from T a p. 
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We now turn to a description of the apparent error in 
the momentum constraint equation used by WMM. Con- 
sider the following two inequivalent definitions of momen- 
tum density. The first is S a = —h^Tp^ri 1 , which is just 
the quantity which appears in the momentum constraint 
Using the perfect fluid stress energy tensor (||) and 



the notations W = 
be written as 



-npu 



P - 



era and a = p + p, it can 



(12) 



The second definition is simply the expression ([12]) with- 
out the projection tensor: 



(13) 



S a = Wau c 



WMM appear to confuse the two different quantities (jl£ 
and (|ll]). They define only a 3- vector Sf, this definition 
[Eq. (47) of Ref. is compatible with both definitions 
(U) and @, since Si = Si (but S t ^ S t ). However, the 
4- vector S a appears in some of their equations. Their 
hydrodynamic equations are correct only if their S a is 
interpreted to be S a , while their momentum constraint 
is correct only if S a is interpreted to be S a . 

This confusion apparently gives rise to an error in their 
equation for the shift vector. WMM solve for the shift 
vector by combining the relation (||) with the assumption 
K = and with the momentum constraint ([|). The 
resulting equation for the shift vector is 



D 2 ^ + -D l (D k p k ) = A ij Dj ln(a/(^ 6 ) + 16natp 4 S 3 , 



(14) 



where Di is the derivative operator associated with the 
flat metric Tq, and A ij = D 1 P 3 + D 3 - 2(D k /3 k )f 3 /3. 
Equation ( |l4|) agrees with WMM's corresponding Eq. 
(33) of Ref. However, WMM then rewrite their vari- 
able S 3 in terms of Sj = Wauj. For the correct variable 
S 3 = S j , we have S j = (p~ 4 Sj = ip- 4 Wau 3 . For the 
incorrect variable S 3 we have instead 



S 3 = Wa 



-4 
(f Uj 



W ■ 
—fi 3 
a 



(15) 



Inserting Eq. (|15|) into Eq. (|1J), WMM obtain the equa- 
tion [Eq. (41) of Ref. @, also Eq. (15) of Ref. |]] 



8tt \^ F 'T a p + AT a p] , where g^ v is the metric obtained 
by solving the WMM equations and ^- > T Q/ 3 is the fluid 
stress-energy tensor (jfy. It is a useful point of view 
to regard Einstein's equation as being satisfied exactly, 
but with an extra type of matter present whose (con- 
served) stress tensor is AT Q/ g and which interacts with the 
neutron stars only gravitationally. The fictitious stress- 
energy tensor is of post-l-Newtonian order, although 
without the error term in Eq. ( |l6| ) it would have been 
of post-2-Newtonian order. Our approach will be to cal- 
culate AT Q( 3 analytically to post-l-Newtonian order, and 
then, starting from a correct, post-l-Newtonian descrip- 
tion of the binary, to solve for the perturbation to the 
stellar structure that is linear in AT a p. 

In our calculations, it will be sufficient to restrict at- 
tention to stationary solutions for which the vector field 
d/dt is a killing vector field, since the numerical, dynamic 
solutions to the WMM equations relax to such stationary 
states [§. There are two contributions to AT a p : (i) a di- 
rect contribution due to the error term in Eq. (|l6|), and 
(ii) an indirect contribution due to the fact that the first 
error causes a non-zero K and invalidates the maximal 
slicing assumption. 

We first calculate the direct contribution. We can de- 
compose the fictitious stress energy tensor as 



AT al3 = Apn a nP 



2n^ a AS f3) + ASg? F) 



-h af3 AS, 



(17) 



where Ap = AT a0 n a n , AS a = -{g afi + n a n /3 )AT f3 



AS 



{g a P + n a nP)AT ap and AS^ F) is orthogonal to 

n a and tracefree. We define A/3 J = (3 l — ^^ftjXk, where 
e^k is the volume element associated with the flat metric 
7ij and x l are Cartesian coordinates associated with 7^ ; 
thus A/3* — > as r — * 00. Rewriting Eq. (|l6|) in terms of 
A/3* and the contravariant components of the 4-velocity 
and taking the post-l-Newtonian limit yields 



D 2 A/3 3 + ^D j (D k Ap k ) = 



(18) 



where pm is the Newtonian mass density and v l — dx l /dt 
is the 3-velocity in the rotating frame (|]) . From Eq. (|l8| ) 
we see that there is there is a direct contribution 



D 2 (3 J + -D 3 (D k p k ) = A 1J Diln(a/tp b ) 



-16ira(p Wa 



W 

Lp Uj - e — (3 3 
a 



(16) 



with £q = 1 ft The correct version of this equation is 
given by e = 0; see j f° r example, Eq. (2.16) of Ref. ||. 

We now turn to calculating the leading order effect 
of this error on the stars' central densities. We define 
a fictitious stress-energy tensor AT a p by G a f3[g^v] = 



e (J~J x x)p M 



(19) 



to the quantity AS 1 . 

Consider now the indirect contribution. Using Eq. (Q) 
and the stationarity condition ip — yields 



K= - iD^ 1 + 6(3 l D l \nL P ] 
a 



(20) 



We now solve Eq. (|T|) for the quantity D^ 1 = D t A(3\ 
insert the result into Eq. (pOl) , make use of the Newtonian 
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continuity equation in the rotating frame Di(pMV l ) = 
—f>M — 0, and use the post-l-Newtonian relation ip = 1 — 
$/2 between the conformal factor ip and the Newtonian 
potential $. The result is 



K = -3e (rj x x) • V$. 



(21) 



Now WMM insert the assumption K = into Eqs. 
(|J), (|^) and (|J). Since K is actually non- vanishing, 
this gives rise to the following contributions to AT a p: 
Ap = K 2 /24tt, AS i = -D'K/Uir, and AS = (K — 
ftDiK - aK 2 /2)/4na. Using the relation @, taking 
the post-l-Newtonian limit, adding the direct contribu- 
tion (|l9|), using the stationarity assumption and letting 
e -> 1 finally yields § Ap = 0, 

AS = —V • [(□ x x) • V$] - (Six x)p M , (22) 

47T 



AS* = — [(O x x) • V] 2 $. 

47T 



(23) 



We next calculate the effect of the fictitious stress en- 
ergy tensor (p2[)-(p3|) on the neutron stars' central densi- 
ties. Focus attention on one of the two stars, say star A. 
We define the two dimensionless parameters e = M/R 
and re = R/L, where M is the mass and R the radius 
of either star, and L is the orbital separation. We will 
work to the leading non-vanishing order in re, which will 
turn out to be linear in re. Now it is known that the 
leading order (tidal) fractional corrections to the inter- 
nal structure of star A due to the other star scale as re 3 , 
to post-l-Newtonian order as well as in Newtonian grav- 
ity 0|; we can neglect these corrections. Hence, accurate 
to 0(re 2 ), we can find a non-rotating coordinate system 
(F, x l ) near star A in which the metric is that of an iso- 
lated neutron star. These coordinates are related to the 
original co-rotating coordinates (t, x 1 ) of the line element 
©by 



t = t + z i ( : £)x i 
R i3 (t)x 3 =z i (t)+x i . 



(24) 
(25) 



Here Rij (i) is the rotation matrix satisfying Ra (t) — 
eijk£ljRki{t), z l (t) = Rij(t)z J A -, and z\ is the (time- 
independent) coordinate location of the center of star 
A in the (t, x l ) coordinates. The transformation (|24|)- 
( p5| ) is approximate but is sufficiently accurate for our 
calculation. 

Next, we combine Eqs. @ and (||)-(j2|) together 
with n M = (1,— P % )/a w (1, — fl x x) to obtain the 
contravariant components of AT"' 3 in the (t, x l ) coor- 
dinate system. We then use the zeroth order metric 
ds 2 = —dt 2 + 2(0 x x)i dtdx 1 + dx 2 to obtain the co- 
variant components AT a p, and finally use the transfor- 
mation (|2^)-(|25|) to calculate the components AT a p in 
the (F, x l ) coordinates, discarding all post-2-Newtonian 



Ap F = - — [(n x x) • V] 2 $ + 2(n x xfpM (26) 

27T 

is the fictitious density, 

A PF = -i- [(O x x) • V] 2 $ + §(n x x) 2 p M (27) 

1Z7T O 

is the fictitious pressure, and where Qjj is a traceless ten- 
sor that does not contribute to the leading order change 
in central density. In deriving Eqs. (26)-(|27|) [but not in 
deriving the expression ( p2|) for ATjj] we replaced SIxza 
by $1 x x, which is valid to leading order in re. 

Consider now the case where star A is non-rotating 
and hence spherically symmetric. Then, the fictitious 
momentum density AT^ will not affect the central den- 
sity of the star [|lj, so we can restrict attention to Ap^ 
and App. To leading order in re, we can replace $ and pu 
in Eqs. (p6|)-(p7|) by the self potential 4>a and the mass 
density pu.A of star A, and we can replace the quantity 
f2 x x by n x za- The first term in Eqs. ( p6| ) is then pro- 
portional to d 2 &A/dy 2 with a suitable choice of y-axis, 
which becomes V 2< &a/3 = 47tpa/,a/3 when we average 
over solid angles about the center of star A. Hence to 
leading order in re we obtain 
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ApF = 3 w orb Pm,a, Ap F = - v 2 lh p M ,A, (28) 

where v or b = O x Z A is the orbital velocity. 

It is clear that the fictitious density and pressure J28j ) 
will cause a fractional increase in the central density of 
star A proportional to ere at leading order. To evalu- 
ate the constant of proportionality we solve for the per- 
turbation to the structure of star A using the following 
modified form of the TOV equations p2|: 



dm 
dr 
dp 
dr 



= 4irr 2 [p + App] , 

(p+p) [m + 4:TTr 3 (p + Ap F )] 
r(r — 2m) 



(29) 
(30) 



terms. The result is AT H = Ap F , AT fl 



-AS 1 and 



ATjj = Appd^ + Q Tv where AS"' is given by Eq. <M 



Our procedure consists of: (i) solving Eqs. (29)-fl30|) 
without the correction terms to obtain the unperturbed 
structure of the star. We use the same stellar model as 
used in Ref. [0, described by the polytropic equation of 
state p = Kp~ M , p = p M + Kp T M /(T - 1), where p M is 
the rest-mass density, T = 2, and K = 1.8 x 10 5 erg cm 3 
gr~ 2 . We choose a central density for the unperturbed 
star of p c = 5.93 x 10 14 gr cm -3 , which implies a bary- 
onic mass of 1.62Mq and a total mass of 1.51M Q . (ii) We 
use Eqs. ( p8| ) to calculate App and App. (iii) We insert 
these App and App into into Eqs. (p9|)-(|3C|), and adjust 
the choice of central density p c until a perturbed stellar 
model with the same total baryonic mass of 1.62M is 
obtained. The result is shown in Fig. ^, where we choose 
the stellar separation to be given by re = R/L = 1/4, 
corresponding to u rb = 0.132. The central density has 
increased by ~ 15%. 
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FIG. 1. Consider a neutron star described by a polytropic 
r = 2 equation of state, at an orbital separation of 4 stellar 
radii from another similar neutron star. The density profile of 
such a star will be very close to that of an isolated star, which 
is shown as the solid line. When the leading order effects of 
the erroneous term in the momentum constraint equation are 
included, the result is the dashed line. 

Turn next to the case when star A is rigidly co-rotating. 
The fractional corrections to its internal structure due to 
its own rotation scale as k 3 , and therefore to leading or- 
der in k the above analysis of the effects of App and 
App is still valid. However, there is now in addition a 
gravitomagnetic interaction between the fluid's velocity 
and the fictitious momentum density AT^. A straight- 
forward computation shows that the radial component of 
the gravitomagnetic force averaged over solid angles is 

-^p M n 2 r[2\(r)+r&(r)}, (31) 

where A(r) is given by (rd r + 3)A = $ with A finite as 
r — > 0. This force gives rise to a fractional change in cen- 
tral density proportional to en 2 . Evaluating this change 
numerically at k = 1/4 for the same stellar model as 
above gives a contribution to 5p c /p c of less than one 
percent. Therefore, the dominant contribution to Sp c 
should be that from App and App, and the crushing ef- 
fect should be seen in the co-rotating case as well as in 
the non-rotating case. 

To conclude, we compare our predictions with the be- 
havior seen in the WMM simulations: (i) The predicted 
magnitude of 5p c / p c agrees with that seen, (ii) The scal- 
ing Sp c oc k oc 1/L is not inconsistent with the scaling 
seen in the simulations Jp| . (iii) Our analysis cannot ex- 
plain the claim by WMM j| that the crushing effect is 
not seen in the co-rotating case. In any case, it should 
be straightforward to verify or falsify our proposed expla- 
nation by re-running the simulations without the extra 
term in the momentum constraint equation. 
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